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A MEETING POINT OF ENTROPY AND BIFURCATIONS IN 
CROSS-DIFFUSION HERDING 

ANSGAR JUNGEL, CHRISTIAN KUEHN, AND LARA TRUSSARDI 


Abstract. A cross-diffusion system modeling the information herding of individuals is 
analyzed in a bounded domain with no-flux boundary conditions. The variables are the 
species’ density and an influence function which modifies the information state of the 
individuals. The cross-diffusion term may stabilize or destabilize the system. Further¬ 
more, it allows for a formal gradient-flow or entropy structure. Exploiting this structure, 
the global-in-time existence of weak solutions and the exponential decay to the constant 
steady state is proved in certain parameter regimes. This approach does not extend to 
all parameters. We investigate local bifurcations from homogeneous steady states analyt¬ 
ically to determine whether this defines the validity boundary. This analysis shows that 
generically there is a gap in the parameter regime between the entropy approach validity 
and the first local bifurcation. Next, we use numerical continuation methods to track the 
bifurcating non-homogeneous steady states globally and to determine non-trivial station¬ 
ary solutions related to herding behaviour. In summary, we find that the main boundaries 
in the parameter regime are given by the first local bifurcation point, the degeneracy of 
the diffusion matrix and a certain entropy decay validity condition. We study several 
parameter limits analytically as well as numerically, with a focus on the role of changing a 
linear damping parameter as well as a parameter controlling the cross-diffusion. We sug¬ 
gest that our paradigm of comparing bifurcation-generated obstructions to the parameter 
validity of global-functional methods could also be of relevance for many other models 
beyond the one studied here. 


1. Introduction 

In this paper we study the following cross-diffusion system: 

(1) dtui = div(VMi - g{ui)Vu 2 ), 

(2) dtU2 = div((5VMi -I- kVu2) + /(mi) — ctU2, 

where Ui = U 2 = U 2 (t,x) for {t,x) G [0,T) x hi, T > 0 is the final time, hi C 

{d > 1) is a bounded domain with sufficiently smooth boundary, V denotes the gradient, 
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div = V- is the divergence and dt = denotes the partial derivative with respect to 
time. The equations are supplemented by no-flux boundary conditions and suitable initial 
conditions 


. . (Vui - g{ui)Vu2) ■ i' 

{6'Vui + kVu2) ■ V 


0 

0 


on t > 0, Mi(0,x) = U 2 { 0 ,x) = in fl, 


where u denotes the outer unit normal vector to dQ. The function Ui{x, t) e [0,1] represents 
the density of individuals with information variable x G at time t > 0, and U 2 {x, t) is an 
influence function which modifies the information state of the individuals and possibly may 
lead to a herding (or aggregation) behaviour of individuals. The influence function acts 
through the term g{ui)'Vu 2 in ([T]). The non-negative bounded function g{ui) is assumed to 
vanish only at Mi = 0 and Mi = 1, which provides the bound 0 < mi < 1 if 0 < mi( 0, x) < 1. 
In particular, we assume that the influence becomes weak if the number of individuals at 
fixed X G is very low or close to the maximal value mi = 1, i.e. 5^(0) = 0 and 5^(1) = 0, 
which may enhance herding behaviour. The influence function is assumed to be modified 
by diffusive effects also due to the random behaviour of the density of the individuals with 
parameter 5 > 0, by the non-negative source term /(mi), relaxation with time with rate 
a > 0, and diffusion with coefficient k > 0. 

If (5 = 0, equations ([I])-([2]) can be interpreted as a nonlinear variant of the chemotaxis 
Patlak-Keller-Segel model |KS70j . where the function U 2 corresponds to the concentration 
of the chemoattractant. The model with nonlinear mobility 5 '(mi) was first analyzed by 
Hillen and Painter |HP02] . even for more general mobilities of the type mi/3(mi)x(m 2). 
Generally, the mobility g{ui) = mi(1—mi) models finite-size exclusion and prevents blow-up 
phenomena |Wrz04] . which are known in the original Keller-Segel model. The convergence 
to equilibrium was shown in |JZ09] . Such models were also employed to describe evolution 
of large human crowds driven by the dynamic field U 2 |BMP11] . 

System 0-0 is one possible model to describe the dynamics of information herding 
in a macroscopic setting. There exist other approaches to model herding behaviour, for 
instance using kinetic equations |DL14j or agent-based models |LS08j . but the focus in 
this paper is to understand the influence of the parameters 6 and a on the solution from 
a mathematical viewpoint, i.e., to investigate the interplay between cross-diffusion and 
damping. 

Equations 0-0 with (5 > 0 can be derived from an interacting “particle” system 
modeled by stochastic differential equations, at least in the case g{ui) = const, (see |GS14] ). 
One expects that this derivation can be extended to the case of non-constant g{ui) but 
we do not discuss this derivation here. The above system with g{ui) = Ui was analyzed 
in |HJllj in the Keller-Segel context. The additional cross diffusion with 5 > 0 in ([2]) was 
motivated by the fact that it prevents the blow up of the solutions in two space dimensions, 
even for large initial densities and for arbitrarily small values of <5 > 0. The motivation 
to introduce this term in our model is different since the nonlinear mobility g{ui) allows 
us to conclude that mi G [0,1], thus preventing blow up without taking into account the 
cross-diffusion term 6Aui. Our aim is to investigate the solutions to ([I])-([2]) for all values 
for (5, thus allowing for destabilizing cross-diffusion parameters 5 < 0. 
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One starting point to investigate the dynamics is to consider the functional structure of 
the equation. In this context entropy methods are a possible tool [JiinlSj . The entropy 
structure can frequently be used to establish the existence of (weak) solutions. Further¬ 
more, it is helpful for a quantitative analysis of the large-time dynamics of solutions for 
certain reaction-diffusion systems; see, e.g., |DF07] . The method quantihes the decay of 
a certain functional with respect to a steady state. An advantage is that the entropy 
approach can work globally, even for initial conditions far away from steady states. More¬ 
over, the entropy structure may be formulated in the variational framework of gradient 
flows which allows one to analyze the geodesic convexity of their solutions [LM131IZM15] . 
However, this global view indicates already that we may not expect that the approach is 
valid for all parameters in general nonlinear systems. Indeed, in many situations, global 
methods only work for a certain range of parameters occurring in the system. The question 
is what happens for parameter values outside the admissible parameter range and near the 
validity boundary. 

One natural conjecture is that upon variation of a single parameter, there exists a single 
critical parameter value associated to a hrst local bifurcation point (5b beyond which a 
global functional approach does not extend. In particular, the homogeneous steady state 
upon which the entropy is built, could lose stability and new solution branches may appear 
in parameter space. Another possibility is that global bifurcation branches in parameter 
space are an obstruction. In our context, the generic situation is different from the two 
natural conjectures. 

In the context of CDi-O. the main distinguished parameter we are interested in is 6 . 
Here we shall state our results on an informal level. Carrying out the existence of weak 
solutions and the global decay to homogeneous steady states 

U = 

via an entropy approach, we hnd the following results: 

(Ml) Using the entropy approach, one may prove the existence of weak solutions to ([I])- 
() 2 |) in certain parameter regimes. 

(M2) The global entropy decay to equilibrium does not extend to arbitrary negative 6. 
Suppose we £x all other parameters, then there exists a critical 6e (to be dehned 
below) such that global decay occurs only for (5>(5e ((5^0). 

(M3) If we consider the limit a —)■ -l-cxo then we can extend the global decay up to 

(5* := —n/'j < 0 , where 7 := max g{v), 

vG[0,1] 

i.e., global exponential decay to a steady state occurs for all (5 > (5* ((5 7 ^ 0) if a; is 
large enough. 

(M4) In the limit o; —)■ 0, we hnd that 6e —)■ -|-cxo. In particular, the entropy method 
breaks down in this limiting regime in the formulation presented here. 

We stress that the results for the global decay (M2)-(M4) may not be sharp, in the sense 
that one could potentially improve the validity boundary 6e- Interestingly, we shall prove 
below that (M3) is indeed sharp for certain steady states, i.e., no improvement is possible 
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in this limit. The proofs of (M1)-(M4) provide a number of technical challenges, which are 
discussed in more detail in Section ETTl and Section [31 We also note that the entropy method 
dehnitely does not extend to any negative It is clear that a global decay to a homogeneous 
steady state for all initial conditions is impossible if bifurcating non-homogeneous steady 
state solutions exist as well. We use analytical local bifurcation theory for the stationary 
problem, based upon a modihcation of Crandall-Rabinowitz theory |Kie04] . to prove the 
following: 


(M5) The bifurcation approach for homogeneous steady states can be carried out as long 
as 

S j^Sd-.= -K/g{ul). 

On a generic open and connected domain, local bifurcations of simple eigenvalues 
occur for 


= <5d + 


/4n - 


Z'K) - 


a 




where /i„ are the eigenvalues of the negative Neumann Laplacian. 

(M6) If a > 0 is sufficiently large and hxed, <5^ < hd < and the bifurcation points 
accumulate at (5d- 

(M7) If a > 0 is sufficiently small and fixed, (5d < bifurcation points again 

accumulate at (5d- 


Although these results are completely consistent with the global decay of the entropy 
functional, they do not yield global information about the bifurcation curves. In general, 
it is not possible to analytically characterize all global bifurcation for arbitrary nonlinear 
systems. Therefore, we consider numerical continuation of the non-homogeneous steady- 
state solution branches (for spatial dimension d = 1). The continuation is carried out using 
AUTO |DCD'*~07] . Our numerical results show the following: 

(M8) We regularize the numerical problem using a small parameter p to avoid higher¬ 
dimensional bifurcation surfaces due to mass conservation. 

(M9) The non-homogeneous steady-state bifurcation branches starting at the local bifur¬ 
cation points extend in parameter space and contain multi-bump solutions, which 
deform into more localized (herding) states upon changing parameters. 

(MIO) A second continuation run considering p —)■ 0 yields non-trivial solutions for the 
original problem. In particular, solutions may have multiple transition layers (re¬ 
spectively concentration regions) and the ones with very few layers occupy the 
largest ranges in h-parameter space. 

Combining all the results we conclude that we have the situations in Figure [11(a)- (b) 
for generic hxed parameter values and a generic hxed domain. These two main cases of 
interest are: 


(Cl) a > 0 sufficiently large: In this limit, the entropy validity boundary, the analytical 
bifurcation approach, and the numerical methods are organized around the singular 
limit at (5 = 5*. Indeed, note that 


6* = 5d, if M = M* maximizes g{u) on [0,1], 
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and we show below that he ^ <5* as a —)■ +cxo. The generic picture for a homoge¬ 
neous steady state so that u\ does not maximize g and a is moderate and hxed is 
given in Figure [U^ a). 

(C2) a > 0 sufficiently small: In this case, the generic picture is shown in Figure [U^b). 
The entropy decay only occurs for very large values h > he- Interestingly, the 
approaches do not seem to collapse onto one singular limit in this case. 

We remark that the condition k, ^ —6g{ui) does not only occur in the numerical con¬ 
tinuation analysis. It occurs in the context of the entropy method as well as the analytical 
bifurcation calculation. It is precisely the condition for the vanishing of the determinant 
of the diffusion matrix that prevents pushing existence and decay techniques based upon 
global functionals further. The condition also prevents analytical bifurcation theory to 
work as the linearized problem does not yield a Fredholm operator. In some sense, this ex¬ 
plains the singular limit as a —?■ -|-cxd from (Cl). Although (Cl) is quite satisfactory from a 
mathematical perspective, one drawback is that the forward problem may not be well-posed 
in a classical sense if 5 < of course, the stationary problem is still well-dehned. 




Figure 1. Sketch of the different bifurcation scenarios; for more detailed 
numerical calculations see Section [5l Only the main parameter 6 is varied, a 
homogeneous branch is shown in black and bifurcation points and branches 
in blue (dots and curves). Only the hrst two nontrivial branches are sketched 
which contain solutions with one transition layer, (a) Case (Cl) with a > 0 
sufficiently large; for a suitable choice of u* and a —)■ -|-oo all three vertical 
dashed red lines collapse onto one line, (b) Case (C2) with a > 0 sufficiently 
small. 


For (C2), we cannot prove sharp global decay via an entropy functional. However, 
the hrst nontrivial branch of locally stable stationary herding solutions can be reached 
in forward time via a classical well-posed problem, and (C1)-(C2) always make sense for 
adiabatic parameter variation. Although we postpone the detailed mathematical study of 
the the limit a —)■ 0 to future work, the observations raise several interesting problems, 
which we discuss in the outlook at the end of this paper. 

In summary, the main contribution of this work is to study the interplay between three 
different techniques available for reaction-diffusion systems with cross-diffusion: entropy 
methods, analytical local bifurcation and numerical global bifurcation theory. Further¬ 
more, for each technique, we have to use, improve, and apply the previously available 
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methods to the herding model problem ([I])-®. Our results lead to clear insight on the 
subdivision of parameter space into regimes, where each method is particularly well-suited 
to describe the system dynamics. We identify two interesting singular limits and provide 
a detailed analysis for the limit of large damping. Furthermore, we compute via numerical 
continuation several solutions that are of interest for applications to herding behaviour 
using a two-parameter homotopy approach to desingularize the mass conservation. From 
an application perspective, we identify herding states with clustering of individuals in one, 
or just a few, distinct regions, as the ones occupying the largest parameter ranges. Hence, 
we expect applications to be governed by homogeneous stationary and relatively simple 
heterogeneous herding states. 

There seem to be very few works |Gabl21 [AAN96j studying the parameter space interplay 
between global entropy-structure methods in comparison to local analytical and global 
numerical bifurcation approaches. Our work seems to be, to the best of our knowledge, the 
hrst analysis combining and comparing all three methods, and also the hrst to consider the 
global-functional and bifurcations interaction problem for cross-diffusion systems. In fact, 
our analysis suggests a general paradigm to improve our understanding of global methods 
for nonlinear spatio-temporal systems, i.e., one major goal is to determine the parameter 
space validity boundaries between different methods. 

The paper is organized as follows. In Section [2l we state our main results and provide 
an overview of the strategy for the proofs respectively the numerical methods employed. 
In particular, the entropy method results are considered in Section 12.11 the analytical local 
bifurcation in Section 12.21 and the numerical global bifurcation results in Section 12.31 The 
following sections contain the full details for the main results. The proofs using the entropy 
method are contained in Section [31 where the weak solution construction is carried out in 
Section 13.11 and the global decay is proved in Section 13.21 Section 0] proves the existence 
of local bifurcation points to non-trivial solutions upon decreasing 5. The details for the 
global numerical continuation results are reported in Section |5l We conclude in Section |6] 
with an outlook, where we discuss several open questions. 

Notation: When operating with vectors we view them as column vectors and use to 
denote the transpose. We use the standard notation for L^-spaces, W^’’^ for the Sobolev space 
with (weak) derivatives up to and including order k in Lf as well as the shorthand notation 

Wk,2 ^ ggg 

|Eva02] for details. Furthermore, ' denotes the associated dual space, when 
applied to a function space. 


2. Main Results 

We describe the main results of this paper, obtained by either the entropy method or 
local analytical and global numerical bifurcation analysis. 

2.1. Entropy Method. First, we show the global existence of weak solutions and their 
large-time decay to equilibrium. We observe that the diffusion matrix of system (fT])-(]2|) 
is neither symmetric nor positive dehnite which complicates the analysis. Local existence 
of (smooth) solutions follows from Amann’s results |Ama89] if the system is parabolic in 
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the sense of Petrovskii, i.e., if the real parts of the eigenvalues of the diffusion matrix are 
positive. A sufficient condition for this statement is <5 > (5d = ~^/ 7 - The challenge here is 
to prove the existence of global (weak) solutions. 

The main challenge of ([I])-(l2]) is that the diffusion matrix of the system is neither sym¬ 
metric nor positive dehnite. The key idea of our analysis, similar as in |HJ11] . is to dehne 
a suitable entropy functional. The entropy is a special Lyapunov functional which pro¬ 
vides suitable gradient estimates. Compared to Lyapunov functional techniques like in 
|Horlll IW 0 IO 2 ] (used for the case 5 = 0), the entropy method provides explicit decay rates 
and, in our case, L°° bounds without the use of a maximum principle. (Note that in 
the system at hand, the L°° bounds can be obtained by the standard maximum principle 
but there are systems where this can be achieved by using the entropy method only; see 
|Junl5j .l For this, we introduce the entropy density 

u2 

h{u) = ho(Mi) + 777-) “ = (“u U2y e [0,1] X M, 

2oo 

where ho is dehned as the second anti-derivative of 1/g, 


(4) 


ho{s) 



s £ ( 0 , 1 ), 


where 0 < m < 1 is a hxed number, and 


5o := 5 if 5 > 0 , So := if — ^/y < 5 < 0 . 

It turns out that the so-called entropy variables w = {wi,W 2 )~^ with wi = h^^Ui) and 
W 2 = U 2 /S 0 make the diffusion matrix positive semi-dehnite for all S>S* := -K/y, 5^0. 
We remark that for 5 = 0 the method does not work and we do not cover this case. In the 
ta-variables, we can formulate ([I])-([ 2 ]) equivalently as 

dfU = div(i?(t(;)Vta) -|- F{u), 

where u = u{w), F{u) = (0, /(mi) — Q!M 2 )''" and 

(^) - (/'r;’) ■ 

The invertibility of the mapping w i-y u{w) is guaranteed by Hypothesis (H3) below. We 
show in Lemma m below that B{w) is positive semi-dehnite if 5 > 5*, 5 7 ^ 0. The global 
existence is based on the fact that the entropy 

(6) H{u{t)) = J^(^ho{ui{t)) + ^^^^^dx 

is bounded on [0, T] for any T > 0; note that we write u = u{t) here to emphasize the time 
dependence of H. A formal computation, which is made rigorous in Section 13.11 shows 
that 
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^0 Jn 


(/(mi) - au2)u2 dx. 


The terms in the hrst bracket dehne a positive dehnite quadratic form if and only if <5 > 5*. 
The second integral is bounded since f{ui) is bounded. This shows that for some ei((5) > 0, 


( 8 ) 


dH 

df 


< -ei(<5) 


Vui 


( 9{ui) 


+ 


|VM2 

“T 


dx + c, 


where the constant c > 0 depends on hi, /, and a. These gradient bounds are essential for 
the existence analysis. 

Before we state the existence theorem, we make our assumptions precise: 

(HI) n C with d^l G a > 0, K > 0, h{u^) G with G (0,1) a.e. 

(H2) / G C°([0,1]) is nonnegative. 

(H3) g G C^([0,1]) is positive on (0,1), ( 7 ( 0 ) = (?(!) = 0, g{u) < 7 for n G [0,1], where 
7 > 0 , and J^ds/g{s) = J^ds/g{s) = +00 for some 0 < m < 1 . 

The condition g{u) < 7 in [0,1] in (H3) implies that (u^ — m)^/(2j) < ho{ui) and hence, 
h{u^) G L^(r2) in (HI) yields G L‘^{Vt) and ^ ^^(H). Hypothesis (H3) ensures that 
the function ho dehned in (jll) is well dehned and of class (needed in Lemma [5]). Its 
derivative h'^ is strictly increasing on ( 0 ,1) with range M, thus being invertible with inverse 
(hg)“^ : M —)■ (0,1). For instance, the function g{s) = s(l — s), s G [0,1], satishes (H3) 
and ho(s) = slogs + (1 — s) log(l — s), where log denotes the natural logarithm. A more 
general class of functions fulhlling (H3) is g{s) = s“(l — s)^ with a, fe > 1. 


Theorem 1 (Global existence). Let assumptions (H1)-(H3) hold and let 5 > —K/ 7 . Then 
there exists a weak solution to ([I])-(l3l) satisfying 0<ui<l in Q,t>0 and 

Ml, U2 G 00; dtui, dtU2 G Tf„^(0, 00; H^{Vt)'). 

The initial datum is satisfied in the sense of . 


We provide a brief overview of the proof. First, we discretize the equations in time using 
the implicit Euler scheme, which keeps the entropy structure. Since we are working in the 
entropy-variable formulation, we need to regularize the equations in order to be able to 
apply the Lax-Milgram lemma for the linearized problem. The existence of solutions to 
the nonlinear problem then follows from the Leray-Schauder theorem, where the uniform 
estimate is a consequence of the entropy inequality (0). This estimate also provides bounds 
uniform in the approximation parameters. A discrete Aubin lemma in the version of |D,T12j 
provides compactness, which allows us to perform the limit of vanishing approximation 
parameters. 

Although the proof is similar to the existence proofs in |HJlll IJiinlS] , the results of these 
papers are not directly applicable since our situation is more general than in [HJlllI.Tiinlh] . 
The main novelties of our existence analysis are the new entropy ([ 6 ]) and the treatment of 
destabilizing cross diffusion (5 < 0). 
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For the analysis of the large-time asymptotics, we introduce the constant steady state 
u* = {ul,U 2 ), where 

u*i = M?, ul = with := — ^ f u^(x) dx, j G {1, 2}, 

a m(i2) 

and m(f2) denotes the Lebesgue measure of f2. Furthermore, we dehne the relative entropy 


H(uju*) = f h{u\u*) dx 

Jq 

with the entropy density 

(9) h{u\u*) = ho{ui\ul) + 7 ^(m 2 - ^^ 2 )^ where ho{ui\ul) = ho{ui) - ho{ul). 

2oo 

Note that ui conserves mass, i.e. Ui{t) := m(f2)“^ da; is constant in time and 

ui(t) = Ui for all f > 0. Thus, by Jensen’s inequality, ho(uilul) > 0. 


Theorem 2 (Exponential decay). Let assumptions (H1)-(H3) hold, let f2 be convex, let f 
be Lipschitz continuous with constant cl > 0, and let 

(10) Jo£i(5) > -clcs, 

a 

where ei(J) > 0 and cs > 0 are defined in Lemmas [J] andlE, respectively. Then, for t > 0, 

(11) H{u{t)\u*) < where x('5) := min ol > 0. 

I Cs CHOq J 

Moreover, it holds for t > 0, 

( 12 ) \\ui{t) - UiWL^in) + \\u2(t) - 

Recall that Jq = if J < 0 and Jq = <5 if J > 0. The values for Jo^i(<^) are illustrated 
in Figure [2l It turns out that ffTOj) is fulhlled if either the additional diffusion J > 0 
is sufficiently large or if y/a is sufficiently small. The latter condition means that the 
influence of the drift term g{ui)Vu 2 is “small” or that the relaxation —au 2 is “strong”. 
The theorem states that in all these cases, the diffusion is sufficiently strong to lead to 
exponential decay to equilibrium. For all parameters hxed, except 6, we conclude from the 
condition (IT 0 |) that there exists a Jg such that exponential decay holds for J > Jg (J 7 ^ 0) 
and we see that 

lim 6e = S* = —n/j 

Q —)-+00 

as a singular limit already discussed above. We remark that the exclusion of the decay for 
J = 0 seems to be purely technical and we conjecture that exponential decay also holds for 
(5 = 0. On the contrary, extensions to (a —)■ 0 are highly nontrivial and we can currently 
not cover this degenerate limiting case using entropy methods. 

Theorem [2] is proved by differentiating the relative entropy H{u\u*) with respect to 
time, similar as in (I^. We wish to estimate the gradient terms from below by a multiple 
of H{u\u*). The convex Sobolev inequality from Lemma [5] shows that the L^-norm of 
g{uiY^'^Vui is estimated from below by /^ho(ni|M*) dx, up to a factor. The L^-norm of 
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Figure 2. Illustration of for k = 1 and 5 = \ (black curves). The 

corresponding singular limit 5* = —k /7 = —4 is also marked (grey dashed 
vertical line). 

Vu 2 is estimated from below by a multiple of f^(u 2 —U 2 )^ dx, using the Poincare inequality. 
However, the variable U 2 generally does not conserve mass and in particular, U 2 ^ U 2 - We 
exploit instead the relaxation term in ([2]) to achieve the estimate 

H{u{t)\u*) + x{S) [ H{u{s)\u*) ds < 0. 

Jo 

Then Gronwall’s lemma gives the result. The difficulty is the estimate of the source term 
f{ui). This problem is overcome by controlling the expression involving f{ui) by taking 
into account the contribution coming from the convex Sobolev inequality. However, we 
need that 6 is sufficiently large, i.e., cross diffusion has to dominate reaction. 

The above arguments hold on a formal level only. A second difficulty is to make these 
arguments rigorous since we need the test function hQ{ui) — hg(M^), which is undehned if 
Ml = 0 or Ml = 1 (since ^ 0 ( 0 ) = —00 and ^ 0 ( 1 ) = +00 by Hypothesis (H3)). The idea 
is to perform a transformation of variables in terms of so-called entropy variables which 
ensure that 0 < mi < 1 in a time-discrete setting. Passing from the semi-discrete to the 
continuous case, the variable mi may satisfy 0 < mi < 1 in the limit. 

2.2. Analytical Bifurcation Analysis. As outlined in the introduction, the hrst natural 
conjecture for the failure of the entropy method is to study bifurcations of the homogeneous 
steady states u* = {ul,U 2 ), which solve 

/.O', 0 = div(VMi - 5 ((mi)Vm2), 

0 = div((5VMi -I- kVu 2) + f{ui) — Q!M2, 

with the no-flux boundary conditions ([3]). To study the bifurcations of u* under variation of 
6 we use the right-hand side of flT^ to dehne a bifurcation function and apply the theory of 
Crandall-Rabinowitz [CR711 [Kie04j . The problem is that u* is not an isolated bifurcation 
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branch as a function of 6 since fixing any initial mass yields a different one-dimensional 
family of homogeneous steady states with 

(14) nl = ^J^Mx)dx>0. 

Hence, the standard approach has to be modihed and we follow arguments that can be 
found in |CKWWl^ ISW091IWX13] . It is helpful to introduce some notations hrst. For 
p > d, let 



A’ 

:= {u G : Vn • = 0 on 

(16) 

3^ 

:= LP(D), 


3^0 

:= {«! G L^(r2) : f^Ui(x) dx = O} 


where the space X includes standard Neumann boundary conditions. Due to the Sobolev 
embedding theorem we know that is continuously embedded in for some 

6 G (0,1). If Neumann boundary conditions hold, then our original boundary conditions (E]) 
hold as well. However, the converse is only true if we can invert the diffusion matrix, i.e., 
as long as 6 ^ 6d = In particular, we shall always assume for the local bifurcation 

analysis of homogeneous steady states that 


(16) 




K 


This implies that me may not hnd all possible bifurcations and the single point when the 
diffusion matrix vanishes has to be treated separately; we leave this as a goal for future 
work. 

Next, we dehne the mapping : A” x df x M —> 3^o x 3^ x ® by 


( div(VMi — g{ui)Vu2) 

(5Ami fi)An2 - au2 + /(mi) 

ui(x) da: — m(r2)M^ 

The hrst two terms are the usual bifurcation functions one would naturally dehne, the 
third term is used to isolate the bifurcation branch for the mapping X, i.e., to avoid the 
problem with mass conservation, while the last two terms take into account the boundary 
conditions. We know that there exists a family of homogeneous steady state solutions 



x(u:,u;,d) = o 

for each <5 G M. The goal is to hnd the parameter values hb such that at 5 = (5b a non-trivial 
(or non-homogeneous) branch of steady states is generated at the bifurcation point; see 
also Figure m We are going to check that X is C^-smooth and the Frechet derivative 
with respect to u at a point u = (hi, U 2 ) is given by 

/rr\ fTj\ (IXUi-X\Y{g'{ui){yu2)Ui +g{ui)'sJU2\ 

(18) Mu) := D„^(h, (5) ( + k^U 2 - aU2 + f{ui)U, 

V In 

where {Ui^U 2 )^ G A’xA’ and As : X xX —)■ 3^o x 3^ x M. We already know from Theorem |2] 
that for all (5 > 5e (<5 7 ^ 0), the homogeneous steady state u* is globally stable. Clearly this 
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implies local stability as well and this fact can also be checked by studying the spectrum 
of As{u*). From the structure of the cross-diffusion equations ([I])-® one does expect 
destabilization of the homogeneous state upon decreasing 6. 


Theorem 3. Let u* = ® homogeneous steady state, consider the generic pa¬ 

rameter case with —n ^ dg{u"{) and suppose all eigenvalues fin of the negative Neumann 
Laplacian on hi are simple. Then the following hold: 

(Rl) (5):A’xT —Fredholm operator with index zero; 

(R2) there exists a seguence of bifurcation points d = 6^ such that dim (5[J')]) = 

1, where A/'[-] denotes the nullspace; 

(R3) there exist simple real eigenvalues An(d) of As{u*), which satisfy A„((5") = 0 . 
Furthermore, An((5) crosses the imaginary axis at with non-zero speed, i.e., 
D^uFiu *^ FfAs^], where 7V[-] denotes the range and span[e{^] = 

The results from (R1)-(R3) hold quite generically (i.e., for 5 ^ 5^ and for generic 
domains [Uhl72] l and yield, upon applying a standard result by Crandall-Rabinowitz 
[CRT 11 ICR731 IKieOdj . the existence of branches of non-trivial solutions 

(Mi[s],M 2 [s], 5 [s]) e T X T X M, (mi[ 0 ],M 2 [ 0 ],( 5 [ 0 ]) = {ul,u* 2 , 5 ^), 

where s G [—SO) -Sq] parametrizes the steady-state branch locally for some small Sq > 0 , and 
(mi[s], M2[s], (5[s]) 7 ^ {ul,U 2 ,d^) for s G [— So, 0 ) U ( 0 ,So]. Slightly more precise information 
about the branch can be obtained using the eigenfunction Cb and we refer to Section 0] 
for the details. The main conclusion of the bifurcation theorem is that we know that the 
entropy method cannot show the decay to steady state for all parameter regions. However, 
to track the non-trivial solution branches in parameter space, it is usually not possible to 
compute the global shape of all bifurcation branches analytically. In this case, numerical 
bifurcation analysis is extremely helpful. 


2.3. Numerical Bifurcation Analysis. The results from Section [2TII2.2I do not provide 
a full exploration of the dynamical structure of the solutions for the parameter regime 
6 < 6*. To understand this regime better we study the bifurcations of ffT3|) numerically for 

(19) f{s) = s(l — s), g{s) = s(l — s), s G H = [0, /] C M. 


for some interval length I > 0. Note that this yields a boundary-value problem (BVP) 
involving two second-order ordinary differential equations (ODEs) 


( 20 ) 0 = 

( 21 ) 0 = 

with boundary conditions 

( 22 ) 0 = 


d 

dx 


dui . , dM2 

- 9{ui) 


dx 


dx 


d^ui d'^U2 , . 

0^^ K-r-ir - 0.U2 /(Ml). 


dx^ 


dx^ 


. (o)-ffMo)-)^(o), 

dx dx 


0 = 5^(0) + k^(0), 


dx 


dx 


0 
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(23) 


dwi, , , , ,, dwo, - 

0 = ^(1 1 , 

dx dx 


^dwi , , duo, s 

0 = S^{1) + K-^{1). 


dx dx 

An excellent available tool to study the problem fl20|) - fl2^ is the software AUTO |DCD^07] 
for numerical continuation of BVPs; for other possible options and extensions we refer to 
the discussion in Section |6l AUTO is precisely designed to deal with BVPs for ODEs of the 
form 


(24) 


dz 

— = F{z-p), 
dx 


xe[0,1], G(ta(0), tc(l)) = 0 


t,N 


t>N 


X —)> R"' , G : X and p G 


^ are parameters and 
as a system in the 


where F : 

z = z{x) G is the unknown vector. It is easy to re-write 
form (12^ of four hrst-order ODEs, i.e., we get iV = 4, consider the scaling x = x/l to 
normalize the interval length to one, then drop the tilde for x again, and let 

Pi := (5, P2 := K, P3 := Pa ■= h 

so P = 4 with primary bifurcation parameter 5. For more background on AUTO and on 
numerical continuation we refer to |KOGV071 IKel771 IGov87j . In the setup fl24|) one can 
numerically continue the family of homogeneous solutions 

{u\5) = {u*i,U2,5) 

as a function of 6, i.e., to compute u* = «*(•; 6) for 6 in some specihed parameter interval. 
Although this calculation yields bifurcation points for some 6 values, it is not straight¬ 
forward to use the formulation fl 20 p - fl^T]) to switch onto the non-homogeneous solution 
branches generated at the bifurcation point. The problem is due to the mass conservation 
since 

/K) 


Ml = m(r2) ^ / Ml dx = M*, 


Mn = 


a 


is a solution for every positive initial mass m^. In particular, the branch of solutions is 
not isolated and there exist parametric two-dimensional families of solutions. There are 
multiple ways to deal with this problem; see also Section [HI One possibility is to resolve 
the degeneracy of the problem via a small parameter 0 < p 1 and consider 


(25) 


(26) 


0 = — 


0 = 


dMi 

dx 


.dV 

0 , „ -I- K 


dx 


9{ui) 


dU2 

dx 


p(mi - Ml), 


d\2 


au 2 + f{ui). 


dx^ dx^ 

for a hxed positive parameter mi > 0. In particular, upon setting 

dMi 


:= Ml, 


^2 := M 2 , 


^3 := 


dx 


2(4 := 


dM 2 

dx 


as well as 


Pb ■= Ml, 


P6 ■= P, 


P = 6, 
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we end up with a problem of the form fl2d|) by transforming the two second-order ODEs 
to four hrst-order ODEs and re-labelling parameters. The vector held for the ODE-BVP 
we study numerically is then given by 

/ PaZ?. \ 


(27) F{z-,p) = 


P4Z3 
PaZa 

Pi[- 9 {zi)f{zi) + P 3 g{zi)z 2 + P 2 g'{zi)z 3 ZA + P2P6{zi - p^)]/'Dg 


V 


/ 


PA[-f{zi) + P 3 Z 2 - Pig'{zi)z3ZA - Pl{zi - P5)P6]/Fg 

where Vg := p 2 + Pig{zi) and the detailed choices for the free parameters are discussed 
in Section [5l Observe that the system fl27)) becomes singular if Vg = 0, which is precisely 
the condition 6 7 ^ —K/g{ui) already discovered above. Therefore, we would need also for 
the numerical analysis a re-formulation (or de-singularization) of the problem to deal with 
this singularity and we postpone this problem to future work. As mentioned above, the 
primary bifurcation parameter we are going to be interested in is 5 = pi. The main results 
of the numerical bifurcation analysis, which are presented in full detail in Section 13 are 
the following: 

(Bl) As predicted by the analytical results, we hnd the existence of local bifurcation 
points on the branch of homogeneous steady states in the parameter region with 
(5 < (5d for the case of sufficiently large a and for <5 > (5d for the case of sufficiently 
small a. At each bifurcation point on the homogeneous branch, a simple eigenvalue 
crosses the imaginary axis. 

(B2) The non-trivial (i.e. non-homogeneous) solution branches consist of solutions of 
multiple ’interfaces’ or ’layers’; branches originating further away from (5d contain 
less layers. The branches can acquire sharper layers upon variation of further 
parameters which is important for information herding. 

(B3) At the local bifurcation points, we observe the emergence of two symmetric branches 
of solutions for the case when the nonlinearities are identical quadratic nonlinearities 
of the form s 1 —)■ s(l — s). 

(B4) We also construct non-homogeneous solutions for p = 0 by a homotopy continuation 
step hrst continuing onto the non-trivial branches in 5 and then decreasing p to zero 
in a second continuation step. 

(B5) Furthermore, we also study the shape deformation of non-trivial solutions upon 
variation of k and the domain length 1. The numerical results show that the main 
interesting structures of the problem have already been obtained by just varying 6 
and a. 


3. Entropy Method - Proofs 

3.1. Proof of Theorem [H First, we prove that the new diffusion matrix B{w), dehned 
in (]5]), is positive semi-dehnite if <5 is not too negative. 

Lemma 4. Assume (H3) and 6 > —k/'j, 6^0. Then the matrix B{w) is positive semi- 
definite, and there exists ei((j) > 0 such that for all z = (zi, ^ 2 )"'" ^ w EMf: 

z'^B{w)z > ei{S){g{ui)zl + zl). 
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It holds ei((5) —)■ 0 as (5 \ 0 and S \ —k /7 (see Figure\^. 



Figure 3. Illustration of ei(5) for k = 1 and 5 = | (black curves). The 
corresponding singular limit 6* = —Kj'y = —4 is also marked (grey dashed 
vertical line). 


For later use, we note that the lemma implies that 
(28) Vw : B{w)Vw > 8,(6) , 

where w = {wi,W2) = (hg(Mi), M2 /<^o) are the entropy variables introduced in the introduc¬ 
tion and A \ B = 'Yhij for two matrices A = {Aij), B = {Bij). 


Proof. Let 2; = {zi,Z2)~^ G Mf. Then 

z'^B{w)z = g{ui)zl - ((5o - 5)g{ui)ziZ2 + Sokz^. 


If 5 > 0, then 60 = 6 and the mixed term vanishes, showing the claim for £i(5) = 
min{l,(5fi;}. If —n/j < <5 < 0, we have (5o = K/7. We make the (non-optimal) choice 


^0 — ^o(^) 



> 0 . 


Then £9 < 1 — (1 — 7^/^)^/4, which is equivalent to (k — 7(5)^ < 4(1 — eo)tP. Thus, using 
g{ui) < 7 (see assumption (H3)), 

z^B{w)z = g{ui)zl - - g{ui)ziZ2 + —zl 

\1 J 7 

2 \ t \ f (^ “ 

= eQg[ui)z^ + (1 - eo)g[ui) ) 
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+ 


1 

7 



{k — 7 ( 5 )^ 
47(1 - £0) 



> eog{ui)zl + 


1 

7 



4(1-£o); 




2 

2 - 


In view of the choice of Eq, the bracket on the right-hand side is positive, and the claim 
follows after choosing £i((5) = min{£o(<5), — (6^ —7<5)^/(4(l — ^o(<5)))]/7} > 0 for —n/'y < 

6 < 0 . □ 


The proof of Theorem [T] is based on the solution of a time-discrete and regularized 
problem. 

Step 1: Solution of an approximate problem. Let T > 0, G N, r = T/N, e: > 0, and 
n G N such that n > d/2. Then ^ Let G be given. 

If /c = 1, we dehne = h'{u^). We wish to hnd G such that 


(29) 


1 

r 


{u{w'^) - u{w’^-^)) ■(j)dx+ I : B{w'^)Vw'^ dx 


+ s 


/ ( + 


■ </)] dx = / F{u{w^)) ■ (j) dx 


for all cf) G where /3 G Nq is a multi-index, is the corresponding partial 

derivative, u{w) = {h')~^{w) for ta G M, and we recall that F{u) = (0,/(mi) — au 2 )~^ ■ By 
dehnition of Iiq, we hnd that ui{w) G (0,1), thus avoiding any degeneracy at mi = 0 or 

Ml = 1. 

The existence of a solution to fl2U]) will be shown by a hxed-point argument. In order 
to dehne the hxed-point operator, let y G M^) and rj G [0,1] be given. We solve the 

linear problem 

(30) a{w,(j)) = G{(j)) for all (/) G iL"(r2; M^), 

where 


a{w,(l)) 

G(0) 


V(j) : B{y)Vw dx -I- e 


I ^ D^w ■ D^(j) + w ■ (j) 

\|/3|=n 


- / {u{y) — u{w^ ^)) dx + rj F{u{y)) ■ (f dx. 

Jn Jn 


dx. 


The forms a and G are bounded on Moreover, in view of the positive semi- 

dehniteness of B{y) and the generalized Poincare inequality (see Chap. 11.1.4 in |Tem97] L 
the bilinear form a is coercive: 


a{w, w) > e 


\D^w\‘^ + Itcp j dx > £c||M;||j:/n(Q) for w G 

m=n ' 


By the Lax-Milgram lemma, there exists a unique solution w G M^) L°°(r2; M^) to 

(1^ . This dehnes the hxed-point operator S : M^) x [0,1] M^), S{y, rj) = w. 
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By construction, S{y,0) = 0 for all y G and standard arguments show that 

S is continuous and compact, observing that the embedding ^ 

is compact. It remains to prove a uniform bound for all hxed points of S{-,ri). Let 
w G be such a hxed point. Then w solves (l30l) with y replaced by w. With the 

test function cf) = w, we hnd that 


(31) 


- / {u{w) — u{w^ ^))-tcdx+ / Vw : B{w)'Vw dx 


+ e 




I dx = ?7 / F{u{w)) ■ w dx. 


Since h'^ = 1/g > 0 on (0,1), Hq is convex. Consequently, ho(x) — ho{y) < hQ(x)(x — y) for 
all X, y E [0,1]. Choosing x = u{w) and y = u{w^~^) and using h'^iniw)) = tc, this gives 


- / {u{w) — u{w^ ^)) ■ w dx > - / [h{u{w)) — h{u{w^ ^))) dx. 

Jn 'T JO 

Since ui = ui{w) G (0,1) and / is continuous, there exists Jm = max^gp,!] f{s) and thus, 

/ F{u{w)) -w dx < / (/m - au 2 )u 2 dx < c/, 

Jn Jn 

where c/ > 0 only depends on fM and a. Hence, (l3T|l can be estimated as follows: 


(32) r] / h{u{w)) dx + r 


/ Vw : BiwWw dx + er ( + 

Jn Jn V 


tap 1 dx 


< yrcf + r] I h{u{w^ ^)) dx 


This yields an bound for w uniform in rj (but not uniform in r or e). The Leray- 
Schauder hxed-point theorem shows the existence of a solution w G M^) to (I5U]) 

with y replaced by w and with rj = I, which is a solution to (1291) . 

Step 2: Uniform bounds. Let be a solution to fl29l) . Set w^'^\x,t) = w^{x) and 
t) = u{w^{x)) for X G H and t E {{k — l)r, kr], k = 1,..., N. At time f = 0, we set 
u'‘^'^)(-,0) = h'Q{u^) and We introduce the shift operator 

for t E {{k — l)r, kr], k = 1,..., N. Then solves 


(33) 

1 

T 


■ (j) dx dt + 


V 0 : dx dt 


+ s 


Jo Jn 


dx dt = 


F(m*^’'^) • 0 dx dt 


for piecewise constant functions 0 : (0,T) -E By density, the weak formulation 

also holds for all ^^(0, T; Tr"(H; M^)). 
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We have shown in Step 1 that the solution w = satisfies the entropy estimate 
By fl25]) . we obtain the gradient estimate 


/ Vw'' : B{w^)Ww‘^ da: > ei(h) min{ 7 -\ / (|V<r + |V<r) da:, 

ao Jo. 

since g{u\) < 7. Thus, we obtain from fl3^ the following entropy inequality: 


(34) 


h{u'^) da: + CqT / da: 


,k\2\ 


+ ST 


L (^ 

V |o|_„ 


) da: < c/r + / h{u^ da:. 


k\2 


k-l\ 


where Cq = £i((5) minjy \ (5o }. Adding these inequalities leads to 

k 

h{u^) da: + cqt 


E /(|V«ir + |Vn^ 2 r) da: 
.7=1 


3= 

k 


+ ST 


^2 [ ( ^2 + da: < Cfkr + f h{iP) da:. 

j=l V \R\=n ' 


Since 


h{u^) dx = i ho(Mi) + 


iu^2? 


da: > / {U 2 Y da:. 


2(^0 / 2 ( 5 o Jn 

the above estimate shows the following uniform bounds: 

(35) ||m[ ^||l°o(0,T;L°°(0)) + ll'W^ ^||l°°(0,T;L2(O)) < C, 

(36) \\u^2\\L2(o,T-,m{n)) + Wu''^\\L2(o,T-,m{n)) < c, 

(37) \\L^{0,T;H”{Q)) < C, 

where c > 0 denotes here and in the following a constant which is independent of e or r 
(but possibly depending on T). 

In order to derive a uniform estimate for the discrete time derivative, let 0 G L^(0,T; 
//"■(hi)). Then, setting Qt = fl x (0,T), 


/ “ Wwi''^)0 da: dt 

T Jq 

X ||V0||i2(Qy) + e||rcj’'^||L2(o,T;ja'"(O))||0||L2(O,T;ja'"(O)) 
< C\/i||0||i,2(o,T;ir"(O)) + c||0|U2(O,T;ja'i(O)), 


< (||VaaS")|U2(Q^) + \\9{u2)\\L^iQr)\\^u2^\\LHQ,)) 


(38) 


(^ 2 ^^ — arU'' 2 ’)(j) dx dt 


(7")^ 


< (<5||VnS")|U2(Q,) + «:||VnrilL2(Q,))||V0|U2(Q, 


h) 


+ ^ll'*^l’"^IU2(O,T;ra-(O))||0||L2(o,T;aa-(O)) + {\\f{u^2)\\L‘2(QT) + 0 '\\u2^\\l^(Qj.))\\4>\\l^(Qj.) 
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r — O'rM*'’'^||L2(0,r;(R'"(U))') < C- 


< C\/i||0||2,2(o,T;/I'‘(U)) + c||0||l2('o^T;R'1(U))) 
which shows that 

(39) 

Step 3: The limit (e,r) 0. The uniform estimates (l36|) and fl39|) allow us to apply the 

discrete Aubin lemma in the version of |DJ12j . showing that, up to a subsequence which 
is not relabelled, as (e, r) —?■ 0, 

(40) —)■ u strongly in L^(0, T; L^(f2)) and a.e. in Qr, 

weakly in L^(0, T; iJ^(f2)), 

dtu weakly in L^(0, T; (iJ"(fi))'), 

—)■ 0 strongly in L^(0, T; if"(fi)). 

Because of the L°° bound (|35l) for (ui’^^), we have 

5'(wi), f{ui) weakly* in L°°(0, T; L“(fi)) 

(and even strongly in U‘{Qt) for any p < oc). Thus, we can pass to the limit (e, r) —)■ 0 


r — cFt-u 


m 


m 


to obtain a solution to 


(5iMl,0) df + 


{dtU2, 0) dt + 


(Vui — g{ui)'\/u2)<p dx dt = 0, 


JQ 

rT 


((5Vmi + kVm2)0 dx df = 


(/(ui) — au2)4> dx df 


for all (j) G L‘^{0,T] H^{Q)). In fact, performing the limit e ^ 0 and then r —)■ 0, we see 
from fl38|) that dtU G L^(0,T; {H^{Q))') and hence, the weak formulation also holds for all 
(j) ^ L'^{0,T] It contains the no-flux boundary conditions ([3]). Moreover, the initial 

conditions are satisfied in the sense of (-ff^(fl; M^))'; see Step 3 of the proof of Theorem 2 
in |Junl5] . This finishes the proof. 


3.2. Proof of Theorem [2]. We recall first the following convex Sobolev inequality which 
is used to estimate the gradient terms in the entropy inequality. 

Lemma 5. Let Q (Z (d > 1) be a convex domain and let (p & be a convex function 
such that 1/0" is concave. Then there exists cs > 0 such that for all integrable functions u 
with integrable 0(m) and 0"(m)|VmP, 

where m(f2) denotes the measure ofQ. 

Proof. The lemma is a consequence of Prop. 7.6.1 in |BGL14] after choosing the probability 
measure dp, = da:/m(f2) on fl and the differential operator L = A —x-V, which satisfies the 
curvature condition CD(l,oo) since r 2 (M) = ^(|V^Mp-|- iVup) > = r(M). Another 

proof can be found in [AMTUOH Section 3.4]. □ 
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Step 1: Uniform bound for the norm of u\. The norm of u\ is not conserved but 
we are able to control its norm. For this, let G be a solution to 

and set = ui{w^). We introduce the notation v = m(r2)“^ f^v(x) da; for any integrable 
function v. This implies that = u^. Employing the test function 0 = (1,0) in 0291) . we 
hnd that — eTw\. Solving the recursion gives 

k k 

u\ = u\ — ST ^ w{= u\ — er ^ w{, 
i=i i=i 

and by fl37j) . we conclude that 

\u^i\t) -ul\ < ^ 

where u^i\t) = u\ for t G {{k — l)r, kr]. Consequently, as (e, r) —)■ 0, the convergence fHOj) 
shows that Ui{t) = for f > 0. 

Step 2: Estimate of the relative entropy. We employ the test function 

0 = (ho(Mi) - h'Q{ul), (m2 - ul)/5o) = {w\ - h'Q{ul),W2 - u^/do) 


in (|2^ to obtain 



For the hrst integral, we employ the convexity of /iq: 


(< - <-^)(^o(<) - KK)) > (ho(nj) - ho(nr^)) - K{ul){ul - 


which yields 


(■4 - 4-')(-4 - 4 ) > -{(4 - 44 - («y‘ - 44 ), 


h>-[ ( 4 ( 4 ) - 4(:4 ')) - hhil f (4 - 4 1) di 

Jn T' Jn 


+ ^ ((“2 - ulY -{u^ ^ - ulY) da;. 

By fl28|) . it follows that 


h > £i{d) 


f (rnu.k4^..k,2 , 

9{u\) dl 


dx = e,{5) / U'K)|V<|^ + 


^0^ 


da;. 
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Lemma [5] then shows that 


{ho{u\) - ho{u\)) dx + 


ei{5) 


Jn 


iVn^r dx. 


cs Jq 

The third integral in fjTTl) is estimated by using Young’s inequality: 


h>^ {Hf-+(«&- dx>-- / 


Summarizing these estimates, we infer from fITT]) that 

[ iho{u’l) - ho{u^~^)) dx - [ {u^ - u^~^) dx 


+ 


1 

ei((5)r 
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[u: 


k-l 


u. 


;)2) dx 


cs Jn 


(ho(ni) - ho(Mi)) dx + 


ei{5)T 


|Vn^|2 dx 


- y X ^ ^(“2)^) ^ ~ «“2)(«2 - U2) dx. 

Adding these equations over k and using the notation as in the proof of Theorem [1] for 
u['^\ we obtain 


(42) 


{ho{u^i\t)) - ho{u^)) dx - h'^iul) / - n?) dx 

Jn 

{{u^ 2 \'t) - ulf - {ul - ulf) dx 


1 

25 o Jn 




{uX>f + (5o ^{ulf) dx ds + — 

oo jQ Jn 



dx ds 

(/(“i^^) - «W 2 ^^)(W 2 ^^ - ul) dx ds. 


<^0 Jo Jn 


Step 3: The limit (e, r) —)■ 0. Because of the L°° bound for it follows that, for a 

subsequence, ui'* Ui weakly* in L°°(0, T; L^(r2)) and thus, as (e:,r) —?■ 0, 


I {u^^\t) — u() dx = / {uY’it) — ul) dx —)■ / {ui{t) — ul) dx = 0, 

Jn J n Jn 

since Mi(t) = ul for f > 0, by Step 1. The weak convergence of {Vu^'^) to Vm 2 in 
L^(0, T; L^(r2)) implies that 


,(t) 


lim inf 

r^O 



0 Jn 


dx ds < / / IVM 2 P dx ds. 


0 Jn 


Furthermore, by the strong convergence u^i'^ —)■ Ui in L^(0, T; L^(r2)), up to a subsequence. 


u 


h) 


—)■ Ul a.e. in Qt = kl x (0,T) and ho{ui) a.e. in Qt- Then the L 
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bound of implies that hQ{u^^^) —?■ ho{ui) strongly in L^(0, T; L^(r2)) for any p < oo. 

Furthermore, we know that —)■ U 2 strongly in L^(0, T; see fHOj) . Therefore, the 

limit (e, r) —)■ 0 in fli2|) leads to 


{ho{ui{t)) - ho(M?)) + ^ j {Mt) - - ulf) dx 


cs Jo Jn 



o(mi) - ho{ul)) dx ds + 


£ 1 ( 5 ) 



IVM 2 P dx ds 


0 Jo Jn 


< 


5, 



0 Jo JQ 


(/(wi) — au 2 ){u 2 — U 2 ) dx ds. 


Now, we estimate the right-hand side. Because of f{uX) = ^^^d the Lipschitz con¬ 

tinuity of / with Lipschitz constant cl > 0, we infer that (recall ([9]) for the dehnition of 
ho(Mi|Mi)) 


J {hoiui{t)\ul) dx - ho(Mi(0)|M*)) + ^ J ((“2(i) - u* 2 f - MO) - 112)^) dx 

ho{ui{s)\ul) dx ds 


+ 


£i(5) 



1 

< — 

<^0 Jo Jn 


< 


cs Jo Jn 



a 


{f{ui) - f{ul)){u 2 - ul) dx ds - — 

Oo Jo Jn 



{u 2 — ^ 2 )^ dx ds 


2(5oci 



(/(wi) - f{u\)) dx ds - 


'0 Jn 


a 



(m 2 — 112 )"^ dx ds 


0 Jo Jn 


< 


2a(5( 



(«! — ul)^ dx ds — 


0 Jo Jn 


a 



{U 2 — U 2 )^ dx ds. 


0 Jo Jn 


Since ui = ul, a Taylor expansion and the assumption l/hQ(ni) = g{ui) < 7 give 



ho{ui\ul) dx ds = 


0 Jn 



0 Jn 


(43) 


(ho(ni) — ho{ul) dx ds 

1 , 



. * \ 2 


where is a number between ui and u^. We conclude that 


ho{ui{t)\ul) dx + ^ [ Mt) - u* 2 )‘^ dx 
2oo Jn 


7ci 


cs 



+ 


a 



(^2 — 112 )^ dx ds < / ho(ni( 0 )|n^) dx-|- 


Jo Jn 

1 


0 Jo Jn 


25 , 


0 Jn 


/io(Mi(s)|M^) dx ds 
( 112 ( 0 ) - ^ 2 )^ dx. 
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and recalling the notation h{u\U) = /io(mi|mi) + (112 

f Uf UMTn J , • f 74 

/ n(M(t)|(7) dx -I- mm <-— 

Jn I cs Oido 

Then Gronwall’s lemma implies that 



-u*y/{26o), 
h{u\U) ds < 



h{u{0)\U) dx. 


H{u{t)\U) = [ h{u{t)\U) da; < H{u{0)\U), t > 0, 

Jn 

where x(^) is dehned in flTTD . Finally, taking into account (03]), we estimate 

h{u\U) > - u\f + ^(m 2 - 4)^ 

which shows flT^ and hnishes the proof. 


4. Analytical Bifurcation Analysis - Proofs 


In this section, we are going to prove Theorem [3l The proofs follow closely ideas pre¬ 
sented for similar systems in |CKWWT^ ISW09] IWX13] , which are fundamentally based 
upon an application of results of Crandall and Rabinowitz |CR711 ICR73] : see also |Kie04] 
for a detailed exposition of the these results. Recall that we dehned the spaces A, (V, (Vo 
in f[T3]) and the mapping 

AiAxAxM—)-A’ oxA’xI^ 

in ffT7|) . A hrst step is to investigate the Fredholm and differentiability properties of A. 
Lemma 6. The mapping A satisfies the following properties: 

(LI) T{ufi 6 ) = 0 for all 6 E M. 

(L2) iF{ui,U 2 ,S) = 0 implies that {ui,U 2 ) solves (113]) . 

(L3) T is -smooth with Frechet derivative given by flTS]) . 

(L4) If u{x) = {ui,U 2 ) is a homogeneous state and Sg{ui) 7 ^ —k then DuJ^{ui,U 2 ,S) is 
a Fredholm operator with index zero. 


Proof. For (LI) recall that u* = was the notation for a homogeneous steady state. 

Regarding (L2), observe that the hrst two components of A are just the steady state 
equations flTH]) . Statement (L3) follows from a direct calculation. The problem is to show 
(L4). We follow the argument given in |CKWWT^ IWX13] and consider 

(44) D„A(hi, U2, 5)(17i, U 2 V = U 2 V + B2{Ui, U 2 V, 


where Si : A x A —)■ A’q x A’ x M is dehned by 

^Af/i - div[g'{ui){Vu2)Ui + g{ui)VU2] 
6 AU 1 + kAU 2 - aU 2 + f{ui)Ui 
0 


(45) 


Si 


Pi 

U 2 


and the mapping S 2 : A x A ^ Ao x A x M is given by 

0 

(46) S 2 " ' 


Pi 

P2 


0 


J^Ui(x) dx 
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We observe easily that 02 : x A" ^ 3^o ^ x is linear and compact. We need an 

ellipticity condition and Bi shonld satisfy Agmon’s condition |SW09j . We have ellipticity 
for Bi (in the sense of Petrovskii |Jan98[ ISWn9] ) if 


(47) 


det 


5 


K J 




7 ^ 0 , 


for all ^ ... ,^d) G R'^\{0}. Compnting the determinant this condition just yields 

0 7 ^ (^1 + • ■ ■ + + <^ 5 '('hi)) if and only if — k 7 ^ Sg{ui) 


and ellipticity in the sense of Petrovskii follows. Moreover we need to verify Agmon’s 
condition at a hxed angle 6 G [— 7 r, 7 r). Using |SW09[ Remark 2.5] with 6 = 7 r/ 2 , one 
verihes computing a shifted determinant similar to the previously computed one in Wf\ 
that Agmon’s condition holds for all values of k. In particular, the ellipticity condition gives 
a restriction on the parameters for the bifurcation analysis and not Agmon’s condition. By 
applying |SW09[ Thm. 3.3] we infer that 

: A X A ^ 3^ X 3^ X {0} 

is a Fredholm operator of index zero. Hence A’o x A’ x {0} = Ti{Bi) © lU, where TZ{Bi) is 
the range of Bi and lU is a closed subspace of 3 ^ x 3 ^ x M with dimlU = dimA/'(Hi) < 00 . 
Consequently, since the hrst component of Bi is in A’qj we have 

3^0 X 3^ X M = TZ{Bi) © Ho © span{(0, 0,1)"^} 

where IHq = {(i/i, hfs, e 1H| Hi{x)dx = 0} and W = Wo + span{(l, 0, 0)}. Then 
dim IH = dim Wo + 1. Thus the codimension of TZ{Bi) in A’o x 3^ x M is equal to dim W = 
dimAA(Hi). Hence, Hi : A x A -)■ 3^o X 3^ X M is a Fredholm operator of index zero for 
5g{ui) 7 ^ —fi). Therefore, is a Fredholm operator of index zero as B 2 is a compact 

perturbation. Hence, the result (Rl) in Theorem [3] follows. □ 


It seems difficult to improve the result to include the degenerate cases when k = —5g{uX) 
as this would require to deal with bifurcation problems with non-elliptic operators. The 
next goal is to apply |SW09l Thm. 4.3]. To do so, we need some additional properties 
of A. In particular, in order that bifurcations occur from the homogeneous steady state 
u* = (uj, M 2 ) we need that the implicit function theorem fails. For the following lemma we 
need to be in the case where each eigenvalue gin of the negative Neumann Laplacian on 
H eigenvalue is simple. For the one-dimensional case this always holds, while for generic 
d-dimensional domains the eigenvalues are also simple [Uhl72] . 

Lemma 7. Suppose the eigenvalues of the negative Neumann Laplacian on H C M'’* are 
simple and 5g{uf) 7 ^ —n. Then there exist bifurcation points at S = 61^^ such that the map 
T satisfies the following properties: 

(L5) the null space Af[DuJ-{u*,6f;)] is one-dimensional, i.e., span[e]j] = A/'[D„J^(m*, dj])]; 
(L 6 ) the non-degenerate crossing condition holds, i.e., 

(48) D 5 .^(m*, <5^)e^ i 77[D„^(m*, 5")]. 
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Proof. We start by proving (L5). By fH5l) . the null space of 5) consists of solutions 

for 


(49) 


AU,-g{ul)AU2 

6AUi + kAU 2 - aU2 + f{ul)Ui 

/ Ui{x) dx 


0 , 

0 , 

0 , 


with no-flux conditions on dQ. For any pair u = {ui,U 2 ) G x df, we can expand ui and 
U 2 as a series of mutually orthogonal eigenfunctions of the following system 


(50) 


—Am = fiu in fl, 

= 0 on dil, 


multiplied by constants vectors. Let /i„ > 0 be a simple eigenvalue of fl50l) and is the 
eigenfunction corresponding to fin normalized by dx = 1. Then we dehne 


We obtain 


Ui := / Mi(x)e^„(x) dx, U 2 := / U 2 {x)e^^{x) dx. 


(51) 


e^^Awi dx = -fin / dx = -finUi, 


/ e^„Au2 dx = -fin / U2e^„ dx = -finU2- 

'n Jn 


Now, by multiplying the hrst two equations of fl4^ by and integrating over fl, using 
the boundary condition and fl5ip . we arrive at the following algebraic system for Ui and 
U 2 : 


(52) 


U, - g{ul)U2 = 0, 

{Kfin + a)U2 - if'iul) - 6fin)Ul = 0. 


If (5 > /'{uf)/fin then the system fl52l) has only the zero solution. In this case, we would 
have A/'[D„J^(m*, 5)] = 0 for all 5. In order to have existence of a non-homogeneous solution 
we necessarily require 5 < f'{u\)/fin- In this case the system fl52l) has a non-zero solution 
if and only if 


(53) 




K 

iM) 



a ' 

iM)- 


<5(1 + 


Pn - 


fK) 


a ' 


Taking 5 = (5((, we can rewrite the hrst two equations of (1491) as the system: 


(54) 


Mf/A 1 (-<j«)rk<) {Ui\ , {u,\ 

[aU^J K + 5gg(ti:)V -Z'K) a J\U-J [u-J 


Using (1531) and computing the determinant and the trace of the matrix A we hnd that 
its eigenvalues are Ai = 0 and A 2 = —fin, where > 0 is a single eigenvalue of the 
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problem fl50l) . Let T be the matrix whose columns are the eigenvectors corresponding to 
Ai and A2 respectively: 

« 9{ul) 

7'K) 1 


T = 


We have 


T-^AT = 

Then, by considering the transformation 


0 0 
0 


(55) 


= T 


-1 


Ui 

U 2 


it follows that the hrst two equations of can be uncoupled and we hnd that 

Ap = 0 in n, 

Ag = finQ in 

^ ^ a [ p{x) dx + g{ul) f q{x) dx = 0, 

io in 

Vp ■ V = Vg ■ = 0 on dVL, 

where the genericity condition —k 7^ Sl^g{ul) is used to obtain zero Neumann boundary 
conditions. Recall that pn is a simple eigenvalue of fl50|) with eigenfunction Observe 
that dx = 0, which implies that p = 0 and g = for some constant C are the 

solutions of (156|) . Therefore, it follows that 

(57) {U,,U2)'^ = Ce,M<)AV- 

This shows that = span[e^^(g(M];), 1)"'"] =: span[e|^]. In particular, the 

nullspace is one-dimensional and the result (L5) follows. 

To prove (L6), we argue by contradiction and suppose that fH5]) is not satished. Hence, 
by computing D{u*, , it follows there exists (p, g) such that 


Ap-g«)Ag = png{ul)e_ 


f^n 


in O, 


(58) 


fi)Ag -I- (5b Ap — aq + f'{u\)p = 0 in O, 
/ p{x) dx = 0, 


Vp •i^ = Vg-i^ = 0 on dVL. 

As in the hrst part of the proof, it is helpful to consider a suitable projection and we dehne 
P and Q as 

F := / p(a;)eb(a;) dx, Q '= 

Jq Jn 

Multiplying the hrst two equations (155]) by and integrating over H and using boundary 
conditions one obtains an algebraic system for P and Q given by 

P-g{ul)Q = -giul), 
ifiul) - S^PnjP - + a)Q = 0. 


( 59 ) 
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By the definition of 6^^ the determinant of the matrix of coefficients on the left-hand side of 
the system is zero. This implies that the inhomogeneons linear system has no solution. 
Hence the system (l58ll has no solutions and the result (1^ in (L 6 ) follows. □ 


Note that (L5)-(L6) are just the results (R2)-(R3) claimed in Theorem [3l By apply¬ 
ing |SWn9l Thm. 4.3] we obtain the existence of a non-trivial branch of solutions. There¬ 
fore, the local dynamics of the problem already shows that the entropy method cannot 
provide exponential decay to a distinguished steady state for all parameter values. 


5. Numerical Bifurcation Analysis - Continuation Results 

In Section [2.11 we proved the existence of a weak solution for 6 > 6* = —k/6 as well 
as global convergence to a steady state for <5 > (5e (<5 7 ^ 0); in addition, converges to 
6* = —k /7 as a ^ -|-cxo and (5e converges to -|-cxo as a —)■ 0. In Section 12^ we showed the 
existence of non-trivial solutions for 6 = where 6^ is defined in and in particular 
could be bigger or smaller than (5d = K/g(ui) depending on a. 

The numerical continuation results presented in this section aim to augment and extend 
these results. To simplify the comparison to numerical results, we focus on the case 

«^ = 1, ^(s) = s(l-s), /(s)=s(l-s), 

which yields the condition 5 > = —4 for the validity of the entropy method for a —>■ -|- cxd . 

As already mentioned, the values for (jjj depend on a, so we are going to study a case with 
a sufficiently large fSection [5.2p and the case with a sufficiently small (Section |5A]). Below 
we are going to define the meaning of sufficiently large and sufficiently small. First, we 
want to compare the values that we obtain for with the numerical results. The analytical 
problem did not include the small parameter p and the introduction of this term has the 
effect of shifting the bifurcation points. 


5.1. Comparison between the values of The formula for given in the equa¬ 
tion fl53p does not consider the additional term p. Introducing this term, we get a new 
formula which reads 


(60) 



( k/i + a){p + p) 


, 1 
<5d H- 

p 


f\<] 


Kp + a 

9{ul) 


ap ■ 
g{u\)p.' 


We observe that the formulas fl5^ and fl 6 Up . due to the presence of the term p, will not 
give the same values 5'^ but the two equations correspond if we take p = 0. We fix the 
following parameter values 


(k, a, /, Ml, p) = (1, 0.2, 20, 0.594, 0.05). 

We are interested in computing the values of and to observe how the parameter p shifts 
the bifurcation branches. 

In Table [U we reported the bifurcation points for n E {1, 2,..., 9} computed with 
the two formulas (l53P and fl60P in comparison to the numerical continuation results using 
AUTO. The values detected using AUTO precisely correspond to the values computed with 
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n 

1 

2 

3 

4 

5 

6 

7 

8 

9 

dMl) 

-45.38 

-14.45 

-8.73 

-6.72 

-5.80 

-5.29 

-4.99 

-4.79 

-4.66 

dm} 

-121.89 

-20.81 

-10.50 

-7.51 

-6.24 

-5.58 

-5.19 

-4.94 

-4.77 

AUTO 

-121.89 

-20.81 

-10.50 

-7.51 

-6.24 

-5.58 

-5.19 

-4.94 

-4.77 

Table ; 

. Comparison between the analytical and numerical 

Difurcation va 


ues. The last two rows compare the numerical and analytical solutions with 
0 < 1 . 


the formula (]60|) as expected while the points are shifted in comparison to the values for 

p = 0. 


5.2. Case 1: a sufficiently large. Recall the formula for 61^ given in fl5^ : 


(5 


n 

b 



a ■ 

IM)-' 


We observe that if a > f'{ul)g{ul) then < (5d and the branches will approach the limit 
value (5d for n ^ oo. Since we are using fl60l) . the condition on a is 


a > - 

L p + /i„ J 

and, in the case of an interval we can compute the eigenvalues /i. So, a sufficiently large 
means 


a > - 


(61) 


Figure m shows a continuation calculation for fixed parameters 


/nvrx 2 

\f’{ul)g{ul) - Kpi 

1 1 ) 

p+(!fr J 


{K,a,l,ui,p) = {p2,P3,Pi,P5,P6) = (1,0.2,12,0.594,0.05) 


using 6 as the primary bifurcation parameter. We observe that the condition on a is 
satisfied since the right-hand side of fIM]) is negative for all n G M and a = 0.2. The steady 
state we start the continuation with is given by 


iul,U 2 ) = {ui,f{ui)/a). 

We begin the continuation at 5 = —25 and we find only one bifurcation point when 6 
is decreasing, i.e. for S < —25. This result is expected since 6^ = —121.889 is the value 
corresponding to the first eigenvalue. Moreover, we do not detect any bifurcations for 
6 > —4 = 6*. The interesting results in the bifurcation calculation in Figure 0] occur when 
we increase the primary bifurcation parameter 6. In this case, several branch points are 
detected, in particular the closer we are to the value 5d, the more bifurcation points are 
found. In Figure IU we have shown the first six branch points detected obtained upon 
increasing 6. The point detected at 5 = —20.8116 corresponds to the second non-trivial 
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Figure 4. Continuation calculation for the system (l24l) with parameter val¬ 
ues {K,a,l,Ui,p) = {P 2 ,P 3 ,P 4 :,P 5 ,P 6 ) = (1,0.2,20,0.594,0.05) and primary 
bifurcation parameter 6. (a) Bifurcation diagram in [6, || 2 ;||i 2 )-space showing 
the parameter on the horizontal axis and the solution norm on the vertical 
axis. Some of the detected bifurcation points are marked as circles (ma¬ 
genta). The last branch point (blue circle) is not a true bifurcation point 
but results from the degeneracy 6 = —K/g{ul) =: (5d. At the other branches 
points (magenta, hlled circles) non-homogeneous solution branches (blue, 
cyan, magenta, green...) bifurcate via single eigenvalue crossing. The value 
6* = —k /7 = —4 is marked by a vertical grey dashed line, (b) Solutions 
are plotted for {x,Ui = Ui{x)) at certain points on the non-homogeneous 
branches; the solutions are marked in (a) using crosses. 

bifurcation branch. There are more and more points as we get closer to 5d- The last point 
detected (in blue) is not a bifurcation point but corresponds to the degeneracy at 

K/g{ul) = -1/(0.594(1 - 0.594)) -4.1466. 

The remaining detected branch points in Figure SJare true bifurcation points. This numer¬ 
ical result is in accordance with the analytical results on the existence of bifurcations in 
Theorem [3l In fact, one can carry out the same calculation as in Section 01 At each bifurca¬ 
tion point, a simple eigenvalue crosses the imaginary axis. One can use the branch switching 
algorithm implemented in AUTO to compute the non-homogeneous families of solutions as 
shown for four points in Figure l4l(a). In Figure l4l(b), we show a representative solution 
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Ui = Ui{x) on each of the four solution families. The solutions are non-homogeneous 
steady states and have interface-like behaviour in the spatial variable. Each family has 
a characteristic number of these interfaces. There are families with even more interfaces 
than the one shown in Figure |l](b4), which can be found upon increasing 5 even further; 
we are not interested in these highly oscillatory solutions here. 



Figure 5. Continuation calculation for the system fl2T|) as in Figure H] with 
a focus on the second bifurcation point (hlled circle, magenta). One can show 
that by using two different local branching directions that two different non- 
homogeneous solution branches (red) bifurcate via single eigenvalue crossing 
but the two branches contain solutions with identical L^-norm for the same 
parameter value. This is a result of a symmetry in the problem, (b) Three 
different solutions plotted in {x,Ui = Mi(a:))-space at the parameter value 
5 = —21.8819. The three solutions are marked in (a) using crosses. 

Another observation regarding the continuation run in Figure 0] is reported in more 
detail in Figure |5] with a focus on the second bifurcation point. It is shown that there 
are actually two different branches bifurcating at the same point with families of non- 
homogeneous solutions that are symmetric. In particular, one non-trivial solution branch 
can be transformed into the other by considering n i-A 1 — u; as an illustration we refer to 
three representative numerical solutions on the three branches originating at the second 
bifurcation point as shown in Figure |5](b). 

5.3. Case 2: a sufficiently small. As specihed in (M7) when a < f'{u\)g{u\) then 
> <5d and this means that the branches will approach the limit value (5d from the right. 
As pointed out in 15.11 the condition on a is more complicated since our model contains p. 
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The condition on a becomes 

^i)^K) - 
P + fin -I ’ 

i.e. we must choose an a which satished the inequality for each single fin- We £x 

= (1,0.001,50,0.211325,0.05) 
for the numerical continuation in this section. 


0 < a < 


n 




Figure 6. Continuation calculation for the system (12T1) with parameter 
values {K,a,l,Ui, p) = {P 2 ,P 3 ,P 4 :,P 5 ,P 6 ) = (1,0.0001,50,0.211325,0.05) and 
primary bifurcation parameter 6. (a) Bifurcation diagram in (5, || 2 ;||i 2 )-space 
showing the parameter on the horizontal axis and the solution norm on the 
vertical axis. The detected bifurcation points are marked as circles (ma¬ 
genta). At the three branch points (magenta, hlled circles) non-homogeneous 
solution branches (blue, cyan, magenta) corresponding to 5)^, 5^, bifurcate 
via single eigenvalue crossing. The value 5* = — = —4 is marked by a 

vertical grey dashed line, (b) Solutions are plotted for {x,Ui = Ui{x)) at 
certain points on the non-homogeneous branches; the solutions are marked 
in (a) using crosses. 

With these values the condition on a is given by 0 < a < 0.0033827 which is satished. 

We also hnd that with our choices 

(5d < 5b < <5* < < <5b < < <5b < < 4, n>6, 
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i.e. there are some bifurcation points which are bigger than 6* and some which are smaller 
but all of them are bigger than (fj- We begin the continuation at 5 = 3 and we detect only 
two more branches when we increase <5 at 5 = 43.4851 and 6 = 9.98041 which correspond 
to (5,^ and 61- We focus on the branches for n G {1, 2, 3,4, 5} such that > S*. This case 
is represented in Figure [6l 



Figure 7. Continuation calculation for the system (1241) with parameter 
values (k, a,/, Ml, p) = {P 2 ,P 3 ,P 4 :,P 5 ,P 6 ) = (1,0.0001,50,0.211325,0.05) and 
primary bifurcation parameter 6. (a) Bifurcation diagram in (5, || 2 ;||i;^ 2 )-space 
showing the parameter on the horizontal axis and the solution norm on the 
vertical axis. Some of the detected bifurcation points are marked as circles 
(magenta). The last branch point (blue circle) is not a true bifurcation point 
but results from the degeneracy 6 = —K/g{ul) =: 6d- At the other branch 
points (magenta, hlled circles) non-homogeneous solution branches (green, 
blue, cyan) bifurcate via single eigenvalue crossing. The value 6* = —ft/y = 

—4 is marked by a vertical grey dashed line, (b) Solutions are plotted for 
(x,Ml = Mi(x)) at certain points on the non-homogeneous branches; the 
solutions are marked in (a) using crosses. 

Numerically we observe that all the branches stop when they reach the critical value 5*. 
Next, we consider n > 6 such that 5d < <5^ < 6* as reported in Figure [71 In this case there 
are two critical values: 5* = —4 (dashed line) and (5d = —6 (blue circle). The branches 
detected for a 6 close to 6* have the same direction as the branches detected for 6 > 6*] 
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but starting from a certain n, in this case n = 8, we notice that the branches change the 
direction. Probably this behaviour is due to the fact that the branches cannot cross the 
value 5 = (5d- We do no detect any branch for 5 < (5d- 

In the range between hd and 6* the branches do not seem to overlap. Numerically, 
one observes that the branches get shorter and shorter due to the numerical continuation 
breaking down as the branches approach hd- Looking at the shape of the solutions in the 
different branches we can observe that they have more and more interfaces as we approach 
the limiting value 5d- Moreover, the solutions inside a hxed branch get sharper and sharper 
peaks along the branch (see for example the cyan branch). 

5.4. Continuation in p. The next question is if we can hnd non-homogeneous steady 
states also for the original problem with p = 0. This can be achieved by using a homotopy- 
continuation idea. 

First, we continue the problem in 6 and compute the non-homogeneous solution branches. 
Then we pick a steady state on the non-homogeneous branch and switch to continuation 
in p while keeping 5 hxed. The results of this strategy are shown in Figure [8] (for a = 0.2) 
and in Figure [U] (for a = 0.001). For the hrst three solutions shown in Figure 111(b), this 
strategy works if we start from a very small p. Figure [Hl(c) shows the solution in the branch 
for a p 7 ^ 0: we notice that the solutions for the case p = 0 keep the non-constant prohle 
as for p 7 ^ 0 yielding relevant herding solutions for applications. 

In the case with a sufhciently small, the strategy works better and we indeed hnd non- 
homogeneous steady states for p = 0 as shown in Figure [9](b). Moreover we can also obtain 
herding solutions. We use the starting parameter values 

(k, a, /, Ml, p) = (1, 0.001, 50, 0.211325, 0.05). 

We start from 5 = 10 and the hrst branch we detect is = 9.98041. Once we are in this 
branch, we continue in p for a hxed <5 (in this case 6 = —9). For information herding models, 
solutions which are of particular importance are those with sharp interfaces between the 
endstates, i.e., the solution is near zero and near one in certain regions with sharp interfaces 
in between. These solutions represent a herding ehect in the sense of sharply split opinions. 
More precisely, they indicate for which values of the information variable x we observe a 
herding behaviour, i.e. a concentration of individuals {u ~ 1) at certain values of x. Figure 
[9](b) shows herding in the interval [0, 0.2] U [0.8,1], while only a few individuals adopt the 
information value in [0.3, 0.7]. 

5.5. Solutions and other parameters. In this section we focus on the case with a 
sufficiently small. We are interested in studying, how the solutions change depending on 
the other parameters k and 1. We fix as starting parameters 

{K,a,l,ui,p) = (1,0.001,50,0.211325,0.05) 

and consider the branch 5^. We study the solutions depending on the different parameters. 
In Figure fTOl we show changes along the branch (which bifurcates at 5 = 9.98041). We 
observe that the shape is the same along the branch but the interfaces sharpen as 6 is 
decreased. 
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Figure 8. Continuation calculation for the system (j24l) starting with the 
same basic parameter values as in Figure 0] but with p = 0.001. We stop the 
continuation at the solution points for a certain 5 (as done in Figure |l](a)) 
and change from <5 as a primary continuation parameter to p as a primary 
parameter with the goal to decrease the parameter to p = 0. The values for 
5 are 5 = —16 for the red branch, 5 = —9.4 for the green branch and 5 = —7 
for the blue one. (al)-(a3) Bifurcation diagrams in (p, || 2 ;||i 2 )-space. The 
starting point for the continuation is at the right boundary where p = 0.001 
and then p is decreased. (bl)-(b3) Solutions obtained on the bifurcation 
branches above at the point p = 0 (points are marked with squares in (al)- 
(a3)). (cl)-(c3) Solutions obtained on the bifurcation branches for the initial 
system with p = 0.001. We can observe that also for p = 0 the solutions 
have a non-trivial herding-type prohle. 


In Figure [TT] we show how the solution changes with the length of the domain. We con¬ 
sider / = 20, / = 50 and I = 100. The branch is detected at 5 = —3.28144, 9.98041,43.4851 
respectively. Since we consider the same branch, the shape does not change and length of 
the domain shifts the bifurcation points and just scales the solution. 

When we change the parameter k, the bifurcation points are also simply shifted. We con¬ 
sider K = 1, K = 5 and k, = 10. The branch 6^ is detected at h = 9.98041, —92.2877, —214.999 
respectively. Moreover, for the hrst case the branches approach the value hd from the right. 
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Figure 9. Continuation calculation for the system starting with the 
same parameter value and as in Figured We stop the continuation at 5 = —9 
(as done in Figure [6](a)) and change from <5 as a primary continuation pa¬ 
rameter to p as a primary parameter with the goal to decrease the parameter 
to p = 0. (a) Bifurcation diagram in (p, || 2 ;||i 2 )-space. The starting point 
for the continuation is at the right boundary where p = 0.05 and then p is 
decreased, (b) Solution on the second branch 6^ of non-homogeneous steady 
states at p = 0 (point is marked with squares in (a)). 





Figure 10. Solutions along the branch 6 '^ for the system fIMll with param¬ 
eter values (k, a,/, Ml, p) = (P 2 ,P 3 ,P 4 ,P 5 ,P 6 ) = (1,0.001,50,0.211325,0.05). 

(a) Solution of non-homogeneous steady states at 5 = 8.72901. (b) So¬ 

lution of non-homogeneous steady states at <5 = 5.76477. (c) Solution of 
non-homogeneous steady states at <5 = 1.548. 

while in the other two cases from the left. As for the previous case we consider three dif¬ 
ferent solutions with (almost) the same norm (163.863 for the case (a), 163.872 for (b) and 
163.911 for (c)). 
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Figure 11 . Solutions in the branch 6 ^ for the system fl2^ with parameter 
values {K,a,ui,p) = {P2,P3,P5,P6) = (1,0.001,0.211325,0.05). (a) Solution 
of non-homogeneous steady states at <5 = —3.5154, I = 20. (b) Solution 
of non-homogeneous steady states at 5 = 8.93964, I = 50. (c) Solution of 
non-homogeneous steady states at <5 = 37.9117, I = 100. 




Figure 12. Solutions in the branch 6^ for the system fl241) with parameter 
values {a,l,ui,p) = (P3,P4,P5,P6) = (0.001,50,0.211325,0.05). (a) Solution 
of non-homogeneous steady states at 5 = 8.72901, k = 1. (b) Solution of 
non-homogeneous steady states at 5 = —92.2877, k = 5. (c) Solution of 
non-homogeneous steady states at <5 = —220.578, k = 10. 


In summary, we conclude that k and I do not seem to be the parameters of primary 
importance in our context as we can re-obtain similar solutions and similar bifurcation 
structures for different values of k and I upon varying 6, a as primary parameters. 


6. Outlook 

So far, relatively little attention has been devoted to the study of the parameter space 
interfaces of different mathematical methods. In this contribution, we have analysed as an 
example a cross-diffusion herding model to understand where, and how, the global nonlin¬ 
ear analysis approach via entropy variables is connected to bifurcation analysis techniques 
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from dynamical systems. We have shown that both approaches enconnter similar prob¬ 
lems regarding the degeneracy of the diffnsion matrix and we were able to cover different 
parameter regimes by combining the resnlts of the two methods. 

This paper is only a hrst starting point. Here we shall jnst mention a few ideas for fntnre 
work. 

The next step is to analyse the regime a —?■ 0 and to check whether the limitation in 
(ITT]) on a can be improved, or not. In this regard, one also has to consider in which sense 
the forward problem shonld be interpreted for moderate and small valnes of a and for 
(5 < (5d. Recent work |Liol5] snggests that one shonld not only nse the notion of Petrovskii 
ellipticity for the stationary problem |SW09j bnt also consider it in the parabolic context; 
see the classical snrvey |AV64j . 

The next step is to expand the approach to other examples. In particular, many reaction- 
diffusion systems as well as other classes of PDEs have natural entropies, which can be 
used to study global existence and convergence properties. In the nonlinear case, one 
frequently can also employ approaches from dynamical systems to understand the dynamics 
of the PDE. Using a similar approach as we presented here could be illuminating for other 
examples. For example, it is natural to conjecture that there are examples in applications, 
which exhibit the following characteristics: 

(Zl) There exists one hxed parameter region in which the entropy method yields global 
decay. Upon variation of a single parameter, the validity boundary of the entropy 
method coincides precisely with an isolated local supercritical bifurcation point. 

(Z2) There exists one hxed parameter region in which the entropy method yields global 
decay. Upon variation of a single parameter, the validity boundary of the entropy 
method does not coincide with a local bifurcation point. Instead, the obstruction 
is a global bifurcation branch in parameter space with a fold point precisely at the 
validity boundary. 

In this work, we apparently found a more complicated case as shown in Figure [H How¬ 
ever, it seems plausible that the cases (Z1)-(Z2) should occur even in classical problems 
without cross-diffusion, i.e. react ion-diffusion equations with a diagonal positive-dehnite 
diffusion matrix. Determining whether this is true for several classical examples from 
applications is an interesting open problem. 

Regarding the entropy method [CJM+Ol] IDFOGj . it would be interesting to investigate 
in more detail parametric scenarios for its validity regime. For example, the question arises 
whether it is possible to hnd criteria for the validity range that are computable for entire 
classes of PDEs. The entropy approach relies on upper bounds. Although the bounds 
we present here turn out to be sharp in the sense of global decay dynamics in a suitable 
singular limit, this may not always be easy to achieve as demonstrated by the a ^ 0 case 
discussed above. It would be relevant to estimate a priori, which regime in parameter space 
one fails to cover if certain non-optimal upper bounds are used. As above, carrying this 
out for several examples could already be very illuminating. 

Regarding the analytical and numerical bifurcation analysis, there are multiple strategies 
to deal with the problem of mass conservation, or more generally with higher-dimensional 
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solution manifolds. For example, one may try to compute the entire solution family of 
steady states parametrized by the mass numerically |Hen021 rPSlSj . which yields a numer¬ 
ical continuation problem for higher-dimensional manifolds and not only curves. Further¬ 
more, we have focused on the numerical problem in the one-dimensional setup and com¬ 
puting the two- and three-space dimension cases could be interesting |Kuel41 IUWR14] . 
Regarding analytical generalizations, a possible direction is to view 6* as a singular limit 
and phrase the problem as a perturbation problem |Ni981 IFif731 IAK15] . 
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